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Abstract 

In this paper we address the challenging problem of multiple source localization in Wireless Sensor 
Networks (WSN). We develop an efficient statistical algorithm, based on the novel application of Sequential 
Monte Carlo (SMC) sampler methodology, that is able to deal with an unknown number of sources given 
quantized data obtained at the fusion center from different sensors with imperfect wireless channels. We also 
derive the Posterior Cramer-Rao Bound (PCRB) of the source location estimate. The PCRB is used to analyze 
the accuracy of the proposed SMC sampler algorithm and the impact that quantization has on the accuracy 
of location estimates of the sources. Extensive experiments show that the benefits of the proposed scheme 
in terms of the accuracy of the estimation method that are required for model selection (i.e., the number 
of sources) and the estimation of the source characteristics compared to the classical importance sampling 
method. 

Keywords: Wireless sensor networks, localization, multiple sources, quantized data. Sequential Monte Carlo 
sampler, Bayesian inference 


I. Introduction 

Wireless sensor networks (WSN) are composed of a large numbers of low-cost, low-power, densely dis¬ 
tributed, and possibly heterogeneous sensors. WSN increasingly attract considerable research attention due to 
the large number of applications, such as environmental monitoring [jT|, weather forecasts |[2l-|f5l, surveillance 
il6l . 171 . health care |j8j, structural safety and building monitoring (9J and home automation |5j], flOl . We 
consider WSN which consist of a set of spatially distributed sensors that may have limited resources, such 
as energy and communication bandwidth. These sensors monitor a spatial physical phenomenon containing 
some desired attribute (e.g pressure, temperature, concentrations of substance, sound intensity, radiation levels, 
pollution concentrations, seismic activity etc.) and regularly communicate their observations to a Fusion Centre 
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(FC) in a wireless manner (for example, as in ifTTl — ifTTl '). The FC collects these observations and fuses them 
in order to reconstruct the signal of interest, based on which effective management actions are made iflOl . 

In this paper, we study the source localization problem which is one important problem that arises in WSN, 
see for instance the overviews in fT8l and |fT9ll . 

A. Existing works on source localization from WSN 

A number of works have addressed different aspects of this source localization problem. For instance in a 
distributed sensor localization problem the work of lfl8ll studied the problem of sensor localizations and the 
accuracy of such estimation in ad-hoc WSN based on measurements and statistical model design for WSN 
measurements such as time of arrival, angle of arrival and received signal strength. The observations utilized 
to make this inference typically come from a WSN in which there is typically a large number of inexpensive 
sensors that are densely deployed in a region of interest (ROI). Generally, this makes it possible to accurately 
perform energy based target localization. Signal intensity measurements arc very convenient and economical 
to localize a target, since no additional sensor functionalities and measurement features, such as direction of 
arrival (DOA) or time-delay of arrival (TDOA), are required. 

Such energy-based methods, based on the fact that the intensity (energy) of the signal attenuates as a 
function of distance from the source, have been proposed and developed in I20l - ll26l . More precisely, i2l1l 
developed a least-square method to perform the task of localization for a single source based on the energy 
ratios between sensors. This was then extended under a Maximum likelihood (ML) based framework for 
multiple source localizations in ll22l . In this second work, the proposed method uses acoustic signal energy 
measurements taken at individual sensors to estimate the locations of multiple acoustic sources. By assuming 
that the number of acoustic sources is known in advance, their estimation approach involved a combination 
of a multiresolution search algorithm and the use of the expectation-maximization (EM) algorithm. 

Flowever, in both ITHl . 1221 . analog measurements from sensors are required to estimate the source location. 
For a typical WSN with limited resources (energy and bandwidth), it is important to limit the communication 
with the network. Therefore, it is often desirable that only binary or multiple bit quantized data be transmitted 
from local sensors to the fusion center (processing node). Motivated by such constraints, several papers have 
more recently proposed source localization techniques using only quantized data Il24ll - ll26l . In ||24l , a ML 
based approach has been proposed by using multi-bit (M-bit) sensor measurements transmitted to the fusion 
center. In j26l , the authors have also developed for the same problem an alternative solution based on an 
importance sampler which was utilized to approximate the posterior distribution of the single source given 
the quantized data. However, in both works, perfect communication channels between sensors and the fusion 
center are assumed. Usually, in a target localization scenario, a large number of sensors are deployed in some 
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area where a line-of-sight between sensors and the FC cannot be always guaranteed. In |[25ll . an extension of 
the ML-based approach previously derived in ll24ll has been proposed in order to incorporate the imperfect 
nature of the wireless communication channels. 

B. Contribution 

In this paper, we generalize previous source localization works by proposing a localization algorithm for an 
unknown number of sources given some quantized data obtained at the fusion center from different sensors 
with imperfect wireless channels. The statistical approach we derive is based on the recent and efficient 
sampling framework known as Sequential Monte Carlo Samplers (SMC Samplers) ll27l . |[28l . and is able 
to estimate jointly the unknown number of sources as well as their associated parameters (locations and 
transmitted powers) by providing all the information included in the approximated posterior distribution. In 
addition, we also derive the PCRB which provides a theoretical performance limit for the Bayesian estimator 
of the locations as well as the transmitted powers of the K sources. We demonstrate that the proposed 
framework provides significant improvement over classical importance sampling type methods that have been 
used for a single source context in lf26l and adapted here for multiple sources. 

II. Problem Formulation 

In this section we first present the system model, and then develop the Bayesian framework for jointly 
estimating the unknown number of sources as well as their locations and transmitted powers. 

A. Wireless Sensor Network System Model 

As illustrated in Fig. Q] we are interested in localizing an unknown number of targets in a wireless sensor 
environment where homogeneous and low-cost wireless sensors are utilized. All the sensors report to a fusion 
center which then performs the estimation of the target locations based on local sensor observations. Sensors 
can be deployed in any manner since our approach is capable of handling any kind of deployment as long 
as the location information for each sensor is available at the fusion center. 

Each target is assumed to be a source that follows the power attenuation model. We thus use a signal 
attenuation model to represent the observed power that is emitted by each target l24l . This signal attenuation 
model is based on the fact that an omnidirectional point source emits signals that attenuate at a rate inversely 
proportional to the distance from the source, for instance a traveling wave that may be propagating through 
ground surface or an acoustic pressure wave traveling through free-space medium. 
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x-coordinate [m] 

Fig. 1: Example of two targets in a grid deployed sensor field. 

In this work, as in lf22ll . we will further assume that the intensities of the K sources will be linearly 
superimposed without any interaction between them. The received signal amplitude at the i-th sensor (i = 
1,..., N) is thus given by 

Si = ai + rii, ( 1 ) 


where the measurement noise term, rq, is modeled as an additive white Gaussian noise (AWGN), i.e., n, ~ 
J\f(0, a 2 ) which represents the cumulative effects of sensor background noise and the modeling error of signal 
parameters (the Gaussian assumption is generally admitted since the central limit theorem could be applied 
on a processed signal resulting on the average of the samples received during a time period). The true signal 
amplitude a* from all the targets is defined as li22l : 

( 2 ) 


di — 


k =1 


i.k 


where P*. denotes the k-th source signal power at a reference distance do- The signal decay n is approximately 
2 when the detection distance is less than 1km ll2Tj . Finally djj. corresponds to the distance between the i-th 
sensor and the k -th target: 





Cx,i)^ T {ilk 


(3) 


where (c Xt i,Cy t i) and (:/•/., y/ ; ) are the coordinates of the i-th sensor and the /,:-th target, respectively. In this 
work, we assume that sensor noises as well as wireless links between the sensors and the fusion center are 
independent across sensors, and that a 2 is known (although it is not required for our proposed approach to 
work - this could be indeed embedded in the parameters to be inferred). 
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Center 



Fig. 2: Illustration of the system model. 


As illustrated in Fig. [2] at each sensor, the received signal is quantized before being sent to the fusion 
center. Quantization is done locally at the sensors in order to decrease the communication bandwidth on 
the sensors thereby reducing energy consumption. The data is quantized using an M-bit quantizer (M > 1) 
which takes values from 0 to 2 M — 1 where L = 2 M is the number of quantization levels. The quantizer of 
the i-th sensor transforms its input s, to its output b t through a mapping: E >—>■ {()..... L — 1} such that 


bi = 


0 0 £ ^ ^ 1 , 1 ) 

f A; 1 E Si <C Xi t 2, 

L 1 Xi L —i E Si <C Aq^, 


(4) 


l T 


with A ifi = —oo and Aqj, = +oo. Let Qk = Pi,xi,yi,..., Pk , xk,Vk denote all the K source locations 
and their associated transmitted powers. Under Gaussian assumption of the measurement noise, the probability 
that b{ takes a specific value l € {0,... , L — 1} is: 


p(bi = i\e K ) = q ( Q f\i+i ^ 


a 


a 


(5) 


where Q(-) is the complementary distribution function of the Gaussian distribution defined as: 

f +°° l t 2 
Q(x) = / ~i=e 2 dt. 

Jx V2^ 


( 6 ) 


Finally, the quantized observation are transmitted to the fusion center through an imperfect channel which 


may introduce transmission errors. Let z = 


Zi,...,Z N 


denote the observations collected at the fusion 


center via independent channels from the N sensors. As in ll 1 1 ft . |[25], | |29l . the probability of a received 
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observation Zi taking a specific value j, given the targets’ parameters, 6k, can be written as: 

L—l 

p(zi = j\0 K ) = ^ p(zi = j\bi = m)p(bi = m\9 K ), (7) 

m = 0 

where pj >m := p(zi = j\ b L = m) represents the transition probabilities of the wireless channel, see |QTJ, l25l . 


Since sensor noises and wireless links are assumed to be independent, the likelihood function at the fusion 
center can be written as: 


N 


p{z\e K ) = Y[p(zi\e K ) 


2—1 

N 

n 

2—1 


L—l 


y: p(zi\bi = m)p(bi = m\0 K ) 


,m =0 


( 8 ) 


Concerning the prior information related to the parameters of interest 6 , we use in this work: 

K 


p(Ok) = Y[p(x k ,yk)p(Pk), 


(9) 


k= 1 


where 


p{xk,Vk ) = A Y, p ), 


( 10 ) 


p(Pk ) =IG(a,b), 

with fi p set as the center of the ROI and = diag( a 2 p x ]) is the covariance matrix which is very 
coarse so that its 99% confidence region covers the entire ROI. ZQ(a,b) corresponds to the inverse gamma 
distribution with a and b being the shape and the scale parameter, respectively. Note that the proposed 
inference algorithm does not require the prior distributions to be Gaussian and inverse-gamma and will work 
with other prior distribution choices as required for a given application. 


B. Multiple Source Localization in a Bayesian Framework 


In this work, we are interested in estimating the unknown number of sources as well as their parameters 
(locations and transmitted powers). This problem can therefore be seen as a joint model selection and parameter 
estimation task. We have a collection of K m;iX competing models {M.k}k&{ i,...,/c max } (corresponding in our 
case to the number of sources in the ROI) and one of them generates the observations obtained at the fusion 
center. Associated with each model, there is a vector of parameters 6\. € 0;., where (-)/, denotes the parameter 

space of the model M.^. The objective is to identify the true model as well as to estimate the parameters, 

r i T 


Ok 


Pi i x\ , y\ ,..., 1 7 , Xk , yk 


, associated with this model. 
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Bayesian inference proceeds from a prior distribution over the collection of models, p(A4 k ), a prior 
distribution for the parameters of each model, p{6 k \Ai k ) and the likelihood under each model p(z\O k , Ai k )- 
In order to solve this joint model selection and parameter estimation under this Bayesian framework, we first 
employ a Maximum A-Posterior (MAP) model selection rule and therefore provides a parameter estimate 
under the selected model. Following the model selection step is the inference step of the model parameters 
which can then be deduced from the posterior distribution associated with the model Aik*, i.e. p(0k*\z, Aik * )• 
This procedure is summarised as follows: 

1) Model selection: 


k* = argmaxjpfjA/ffclz)} 
k 

= argmax {p(z\A4k)p(A4 k )} , 

k 

2) Model parameters estimation via Bayes rule: 


( 11 ) 


p(0 k *\z, Ai k *) = 


p{z\0 k *,Ai k *)p(0 k *\A4 k *) 


p(z\M k *) 

Deriving the expressions in (II1H12I) involves calculating the evidence of the A>th model Ai k : 
p(z\M k ) = [ p{z\Qk,Mk)p{Ok\Mk)dO k 


( 12 ) 


'©ic 


N 


n 

i— 1 
N 

n 


L—l 


Y, p{zi\bi = m)p(bi = m\d k ) 
i =1 Lm=0 


1^[ p(x n , y n \Mk)p(Pn\M k )de k 


n =1 


©fc j = i 


L—l 


^ ^ Pi,m f Q 


,m =0 


i,m 


a 


-Q 


^i,m +1 0>i 


a 


nv 

n— 1 


Xn 

Vn 


! Mp! I 




ZQ(P n -,a,b)dG k 

(13) 


Unfortunately, owing to the highly nonlinear observation function of the parameters of interest in Equations 
(HE), the integral in (1 1 31 ) is intractable. As a result, Mk € {1,..., I\ m ax }, both the evidence p(z\Ai k ) and the 
posterior distribution of the parameters p(6 k \z, Ai k ) are intractable. In this work, we propose to use SMC 
sampler in order to have an accurate approximation of both quantities. 


III. Proposed Bayesian Algorithm to Multiple Source Localization in WSN 

In this section we first introduce the general principle of SMC samplers, then develop the SMC sampler 
for multiple source localisation, and finally we derive the point estimate for the parameters of interest. 
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A. General Principle of SMC Samplers 

Sequential Monte Carlo (SMC) methods are a class of sampling algorithms which combine importance 
sampling and resampling. They have been primarily used as “particle filters” to solve optimal filtering 
problems; see, for example, lf30l and OTt for recent reviews. In this context, SMC methods/particle filters 
have enjoyed wide-spread use in various applications (tracking, computer vision, digital communications) 
due to the fact that they provide a simple way of approximating complex filtering distribution sequentially 
in time. However in l27l . |[28ll . there have been a range of developments to create a general framework that 
allows SMC to be used to simulate from a single and static target distribution, thus becoming an interesting 
alternative to standard MCMC methods as well as to population-based MCMC algorithms. Finally, let us 
note that there exists a few other SMC methods appropriate for static inference such as annealed importance 
sampling ||32l , the sequential particle filter of |33l and population Monte Carlo 041 but all of these methods 
can be regarded as a special case of the SMC sampler framework. 

The SMC sampler is based on two main ideas: 

a) Rather than sampling directly the complex distribution of interest, a sequence of intermediate target 
distributions, {77}^, is designed, that transitions smoothly from a simpler distribution to the one of 
interest. In Bayesian inference problems, the target distribution is the posterior 7 tt(0) = p(0\z), thus a 
natural choice for such a sequence of intermediate distributions is to select the following |32|| 

M{0) = oc p(0)p( z|0)^ f (14) 

A 

where {ft} is a non-decreasing temperature schedule with 7 o = 0 and fr = 1 and 77 (0) corresponds 
to the unnormalized target distribution (i.e. 7/(0) = p(0)p(z\0 f ), 'j and Z t = f & p(0)p(z\0)^ >t d0 is the 
normalization constant. We initially target the prior distribution 7r 0 = p(0) which is generally easy to 
sample directly from and then introduce the effect of the likelihood gradually in order to obtain at the 
end, t = T, the complex posterior distribution of interest ttt(0) = p(Q |z) as target distribution. 

b) The idea is to transform this problem in the standard SMC filtering framework, where the sequence of 
target distributions on the path-space, denoted by {77 } f = i , which admits 77(.17) as marginals, is defined 
on the product space, i.e., supp(7r f ) = 0 x © x ... X0 = 0*. This novel sequence of joint target 
distributions 77 is defined as follows: 

Srt(0i -.t) = ( 15) 

where 

t-i 

7t( 0 l:i) = n C k (0 k+1 , 0fc)> (16) 

k= 1 
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in which the artificial kernels introduced {jCfcjjjZ^ are called backward Markov kernels since Ct(0t+i ■ @t) 
denotes the probability density of moving back from Qt+i to 0 t . By using such a sequence of extended 
target distributions {Ttt}J=i based on the introduction of backward kernels {Ci-\ l k J ] , sequential impor¬ 
tance sampling can thus be utilized in the same manner as standard SMC filtering algorithms. 

Within this framework, one may then work with the constructed sequence of distributions, 7 ft, under the 
standard SMC algorithm 1351 . In summary, the SMC sampler algorithm therefore involves three stages: 

1) Mutation: , where the particles are moved from 0 t - 1 to 0 t via a mutation kernel KLt{0t- 1 , 0t) al so called 
forward kernel ; 

2) Correction: , where the particles are reweighted with respect to 7iy via the incremental importance weight 
(Equation (l20l) l: and 

3) Selection: , where according to some measure of particle diversity, such as effective sample size, the 

weighted particles may be resampled in order to reduce the variability of the importance weights. 

In more detail, suppose that at time t — 1, we have a set of weighted particles \ 0^}_ 1 ,W^\ that 

1 J m= 1 

approximates 7r t _i via the empirical measure 


N 


= E W i-i 


(17) 


m= 1 


These particles are first propagated to the next distribution jr t using a Markov kernel K, t {Q t -i,Qt) to obtain 

f (m) 1 N 

the set of particles < 0\ .f > . Importance Sampling (IS) is then used to correct for the discrepancy between 

1 ' Jm= 1 

the sampling distribution rjt(0i : t) defined as 




Hi 


j (m) 

t- 1) °t ) 


( 18 ) 


k =2 


and TTt(0i:t)- In this case the new expression for the unnormalized importance weights is given by 


W, 


( m) 


OC 


<*(» S?) 


tttwt ' 1 1 s —— 1 




(m) a (m) u ,r(m) 


»i(e! m) )nL 2 W(er,.e t , 

where wt , termed the (unnormalized) incremental weights , are calculated as 

kn( rn )\r /n( m ) 

w t (0pl,ol m) ) = ^ t ’ l - l) 


(19) 


( 20 ) 


.. taC n )\r aHi' 

However, as in the particle filter, since the discrepancy between the target distribution if t and the proposal 
r) 1 increases with t, the variance of the unnormalized importance weights tends therefore to increase as well, 
leading to a degeneracy of the particle approximation. A common criterion used in practice to check this 
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problem is the effective sample size 


't — 


N 


Ew 

m= 1 


which can be computed by: 

r > ( £ 

\m=l / 


( m ) \2 




( 21 ) 


If the degeneracy is too high, i.e., the ESS/ is below a prespecified threshold, ESS, then a resampling step is 
performed. The particles with low weights are discarded whereas particles with high weights are duplicated. 
After resampling, the particles are equally weighted. To sum up the algorithm proceeds as shown in Algorithm 

ffl 


Algorithm 1 Generic SMC Sampler Algorithm 
1: Initialize particle system 

2: Sample j f )\ rn ' 1 1 ~ Wi{') an d compute wj m> 

l J m= 1 

ESS 



V-Al 7l(®i 0) ) 

VuiieH) J 



and do resampling if ESS < 


3: for t = 2,. .., T do 

4: Mutation: for each m = 1,..., N : Sample 6™ ~ ; •) where /Ct(s ■) is a 7r t (-) invariant Markov kernel. 

5: Computation of the weights: for each m = 1 ,,N 


W t (m) = 




Normalization of the weights : wj m> = Wl™' 1 



-l 


6: Selection: if ESS* < ESS then Resample 

7: end for 


Let us mention two interesting estimates from SMC samplers. First, since jx t admits tt/ as marginals by 
construction, for any 1 < t < T , the SMC sampler provides an estimate of this distribution 


N 


tt ?{dO) = Y J Wl m) Se< r AdO) 


( 22 ) 


m=l 


and an estimate of any expectations of some integrable function </?(■) with respect to this distribution by 


N 


v = E w t^ e 


t ) 


(23) 


m=l 


Z t f j t (0)d0 

Secondly, the estimated ratio of normalizing constants -= —4 -is given by 


Z t - 


t -1 


/ 7t_i(0)d0 


Z t 

Z t -1 


N 


E wAAmEA” 1 ) 


(24) 


m= 1 
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z t 

Consequently, the estimate of — is 


Zt_ 

Zi 


n 


N 


IIE w£1moE\a 


z, 


k -1 


M ni m )\ 

-l’“k J 


(25) 


k =2 fc=2m=l 

If the resampling scheme used is unbiased, then (l25l) is also unbiased whatever the number of particles used 
m • Moreover, the complexity of this algorithm is O(N) and it can be easily parallelized. 

To conclude this section, let us summarize the advantages of SMC samplers over traditional and population- 
based MCMC methods. First, unlike MCMC, SMC methods do not require any burn-in period and do not 
face the sometimes contentious issue of diagnosing convergence of a Markov chain. Secondly, as discussed 
in 071 . compared to population-based MCMC methods, the SMC sampler is a richer class of method since 
there is substantially more freedom in specifying the mutation kernels in SMC: kernels do not need to be 
reversible or even Markov (and hence time adaptive). Finally, unlike MCMC, SMC samplers provide an 
unbiased estimate of the normalizing constant of the posterior distribution which will be one of the quantities 
of interest in the inference problem tackled in this paper related to finding the number of targets that are 
present in the ROI. 


B. Proposed SMC Samplers for Bayesian Multiple Source Localization 

Since the evidence of the model M. k corresponds to the normalizing constant of the posterior distribution 
of the parameters associated with this model, i.e.: 

ini xz \ _ p{z\0 k ,M k )p(0 k \Mk) _ p(z\e k ,M k )p{e k \M k ) ^ 

fcl p(z]Mk) f ekP (z\e k ,M k ) P (e k \M k )d6 k ’ 

we propose to use the following procedure: 

1) For each model M k , k € 1,... ,K max : approximate the conditional parameter posterior distribution 
p(6 k \z, M. k ) as well as the marginal likelihood p(z\M. k ) using /v rnax SMC sampler algorithms in 
parallel (one for each model {-M-k}kG{i , ax }) using Equations (1221) and (|25T ). respectively. 

2) Perform the MAP Model selection rule: 

k* = argmax{p(z\M k )p{M k )} , (27) 

k 

with p(z\M. k ) corresponds to the unbiased estimate obtained from the A'-th SMC sampler. 

3) Provide a parameter estimate under the selected model, e.g. the minimum mean square (MMSE) 
estimate, using the empirical approximation of the posterior distribution p(0 k *\z, M. k ,) given by the 
fc*-th SMC sampler. 


April 23, 2015 


DRAFT 





12 


Let us now describe in more details the different choices required in the design of each SMC sampler: the 
appropriate sequence of distributions {'ttt}l<t<Ti the choice of both the mutation kernel {K,t}2<t<T and the 
backward mutation kernel {C t -i}t= 2 - 

1) Sequence of distributions 7 ^: An annealing procedure which progressively introduces the effect of 
the likelihood function is chosen as the sequence of intermediate target distributions, i.e. for the k-th SMC 
sampler dealing with model Aik- 



(28) 


where {fk.t} is a non-decreasing temperature schedule with fk.o = 0 and 0k,T = 1- The question that arises 
is how to choose this non-decreasing temperature schedule {fk,t}t=i,...,T- Several statistical approaches have 
been proposed in order to automatically obtain such a schedule via the optimization of some criteria, which 


are known as on-line schemes. fl38l proposed an adaptive selection method based on controlling the rate of the 


effective sample size (ESS^ t), defined in (l2Tb . This scheme thus provides an automatic method to obtain the 
tempering schedule such that the ESS decays in a regular predefined way. However, one major drawback of 
such an approach is that the ESS/,-./, of the current sample weights corresponds to some empirical measure of 
the accumulated discrepancy between the proposal and the target distribution since the last resampling time. 
As a consequence, it does not really represent the dissimilarity between each pair of successive distributions 
unless resampling is conducted after every iteration. 

In order to handle this problem, ll39l proposed a slight modification of the ESS, named the conditional 
ESS (CESS), by considering how good an importance sampling proposal itk,t-i would be for the estimation 
of expectation under Hkj- At the t-th iteration, this quantity is defined as follows: 



N 


CESS M = J2 NW k}-i 


(29) 


In this work, this CESS proposed in |39ll will be used in all the /v max SMC samplers that are run in parallel 
for each model in order to have an automatic specification of their individual temperature scheduling process. 
Owing to the on-line nature of this CESS-based strategy, the total number of iterations performed by each 
sampler is not fixed and does not required to be specified prior to the simulation. 

2) Sequence of mutation kernels Ki-kj: The performance of SMC sampler depends heavily upon the se¬ 
lection of the transition kernels {Ak,t}J = o and the auxiliary backward kernels {£k,t-i}J =2 - There are many 
possible choices for lCk } t which have been discussed in ll27l . Il28l . In this study, we propose to employ MCMC 
kernels of invariant distribution tt;,../ for /Q. ( . This is an attractive strategy since we can use the vast literature 
on the design of efficient MCMC algorithms to build effective and efficient importance distributions fl40l . 
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More precisely, in this work, since we are interested in a complex model with potentially high-dimensional 
and multimodal posterior distribution, a series Metropolis-within-Gibbs kernels with local moves fiOl will be 
employed in order to successively move the B k sub-blocks of the state of interest, Q k) t = Qk,t, 2 , ■ ■ ■ , Qk,t,B k \- 

r’ -| T 

In this work, the sub-block corresponds to the parameters associated to one target, i.e. Qk,t,b = I J i,. .iy, y\ t 
for b = 1,..., £>/,. = k. A random walk proposal distribution is used for each sub-block with a multivariate 
Gaussian distribution as proposal at the i-th iteration of the forward kernel: 


Bk,t,b — Qk,t,b + £ bi 


(30) 


in which e' l b is a Gaussian random variable with zero mean and 3x3 covariance matrix X. The Metropolis 
within Gibbs used in the implementation of the SMC sampler in this paper is summarized in Algorithm [2j 


Algorithm 2 Metropolis-within-Gibbs Kernel •) for the m-th particle 


6 

7 

8 
9 

10 

11 

12 

13 

14 


Initialization Set [g° k t l ,..., el, t ,k\ = #M-i 
for i = 1,..., A'mcmc do 
for b = 1,..., k do 

Sample Q * k t b ~ M s) 

Compute the Acceptance ratio: 


„ (f) . ^ = min p(zK,M k )^<p(e* k ) \ 

( fc ’ fc) ^ h P (z\e k ,M k )^p(0 k )i 


with 


SI 


> 8k,t,b- 1 ’ Qk,t,b+ 1 ’ ■ 

r i i z—1 i —1 i —1 n 

‘ ? ^6 ’ £fc,t,6+l> * ' ‘ fcJ 

Sample random variate u from U(0, 1) 

if u < a(0 k ,d k ) then 

8k,t,b = 8k,t,b 

else 

0k,t,b = 8 l k,t,b 

end if 
end for 
end for 


2 — 1 ] 
> &k,t,k\ 


Set the new particle value at time t as 0 k t = 


„-/Vmcmc ] 

' Uk,t,k J 


and 


3) Sequence of backward kernels £ k ,t-' The backward kernel C k ,t is arbitrary, however as discussed in 
lf28l , it should be optimized with respect to mutation kernel JC k} t to obtain good performance. In |[27| l, |[28|, 
it was established that the backward kernel which minimizes the variance of the unnormalized importance 
weights, W k t , are given by 


a \ _ Vk,t{^k,t)^k,t+l{^k,tAk,t+l) 

"k,t) — ~ 7/1 7 • 


(31) 


Vk,t +1 ( @k,t+ 1 J 

However, it is typically impossible to use these optimal kernels as they rely on marginal distributions rik,.t(@k.t) 
which do not admit any closed form expression, especially if an MCMC kernel is used as t which has 
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a /-invariant distribution. Thus we can either choose to approximate or choose kernels Cj-,t so that 
the importance weights are easily calculated or have a familiar form. If an MCMC kernel is used as forward 
mutation kernel, the following C/,.j is employed 




ttk,t{0k,t—l )%'k,t(0k,t—1 ) @k,t) 


(32) 


7T k,t i@k,t) 

which is a good approximation of the optimal backward kernel if the discrepancy between and 
is small; note that ( 1321) is the reversal Markov kernel associated with /CIn this case, the unnormalized 
incremental weights becomes for the SMC sampler associated to the k’-th model becomes 

roim ) \ 


a( m )\ — 




P (z\e^_ v M k )^-^ 




(33) 


where p(z\0^_ 11 Aik) is defined in Eq. (HI). Expression (l33l) is remarkably easy to compute and valid 
regardless of the MCMC kernel adopted. Note that fk,t. — 4>k,t -1 is the step length of the cooling schedule of 
the likelihood at time t for the k-th sampler. As this step becomes larger, the discrepancy between tt/.^ and 
7 Tfc t_i increases, leading to an increase in the variance of the importance sampling approximation. Thus, it is 
important to construct a smooth sequence of distributions {^k,t} 0<t<T by judicious choice of an associated 
real sequence {4>k,t}J = o as discussed in the previous section. 

Let us remark that when such backward kernel is used, the unnormalized incremental weights in (l33l) at 
time t does not depend on the particle value at time t but just on the previous particle set. It is known that 
in such cases, the particles j 0^ | should be sampled after the weights have been computed and 

after the particle approximation j has possibly been resampled. 


Based on this discussion regarding the different choices, the SMC sampler used in this paper is summarized 
in Algorithm [3] 


C. Point Estimate for the parameters of interest 


Once the model has been selected using the MAP criterion described in (|27T ). the MMSE estimate of the 
parameters for the k *-th model is obtained using ( |23T ): 


N 


Oh* = 


[*>(»)] = E» 


T/J 


m— 1 


(34) 
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Algorithm 3 k- th SMC Sampler Algorithm targeting the posterior distribution p(6 k \z, M k ) 

1: Initialize particle system (set t = 1 and <p k .\ = 0) 

2: Sample { 0^ ) ~ p{0k\Mk) and set = 1/N 

3: while <pk,t < 1 do 

4: t = t + 1 ^ _ 

5: Automatic Cooling procedure: Use a binary search to find cf>* such that CESS^ t = CESS and set cj)k,t = </>* if 

(f>* < 1 otherwise <j)k,t = 1 

6: Computation of the weights: for each m = 1,..., N 


Wt7 = 


Normalization of the weights : 


r {m) 




n -1 


rU) 

k,t 

7: Selection: if ESSk,t < ESS then Resample 

8 : Mutation: for each m = 1,..., N : Sample 0^" 1 ~ •) where ICk,t(-', •) is a 7 Tk,t(-) invariant Markov 

kernel described in more details in Algo. U 

9: end while 

10 : Output: 


11 : Unbiased approximation of the marginal likelihood 

n —2 m =1 


p{z\M k ) 


A %k,n 

11 ry 

n —2 ^k.n—l 


N 


12 : Empirical approximation of the posterior distribution p(9 k \z, Aik) ~ 7r j^ t {d0k) = 

TO= l ’ fc ’* 


where T denotes the last iteration of the &:*-th SMC sampler, since in this last iteration, the system of weighted 
particles represents an empirical approximation of the target posterior distribution, i.e.: 

N 

p(O k ,\z,M k .) = 7r fc . iT (0 fc .) = £ W™6 o g) T {d0 k .). (35) 

m= 1 

Unfortunately, owing to the non-identifiability of the target label in the likelihood and to the same prior 
for each target, the posterior distribution will be multimodal (as it will be illustrated in Fig [5]). The posterior 
is indeed invariant under the permutations of source parameters, i.e., 

p{d k *\z,M k *) =p{ti(0 k *)\z,Mk *) (36) 

where $(■) € V denotes any the permutation for which the posterior is invariant and V is the set of these 
permutations. 

In such a case, the MMSE estimate would lead to very poor performance if selected as a point estimate of 
the source parameters. The problem of having a Monte-Carlo algorithm that approximates such a multimodal 
target posterior, which is invariant under permutation, is known in the literature as the label switching problem 

m. 
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There exists many algorithms that have been proposed in order to deal with this label switching problem in 
the class of Monte-Carlo algorithms. A recent and detailed review of these techniques can be found in li42l . 
Here, we are interested in only post-processing technique in order to extract an accurate point-estimate of 
the state of interest from our particle approximation of the posterior distribution. One of the most commonly 
used relabeling algorithms is the one proposed in flTft . 

Let us denote the unweighted set of particles obtained at the last iteration of the SMC sampler that targets 
the posterior distribution of the selected model by Q = j 0 ^),..., 0'^ j. In the algorithm proposed in BTl . 
one performs inference tasks (e.g. point estimation) as usual but with the relabeled samples, defined as: 

0(fi) = (r?i(0 {1) ),...,i?iv(0 {JV) )), (37) 

where 

■d = (i?i,... ,#at) = arg max L(f?, ■&), (38) 

Vx-xV 

and L(-) is a user-defined cost-function, which is generally chosen as: 

N 

L(B,0) = n^(^(0 W )|/4,Eiv) , (39) 

2=1 

with 

1 N 

A = vE*( 8( '’)- < 40 > 

2=1 

1 N 

^ (41) 

2=1 

The Gaussian cost function in (l39l ) imposes the idea that one wants a relabeled sample to be the most Gaussian 
possible among its permutations $(<2), 0 G V N , in order for 'd(Q) to look as unimodal as possible. 

However, this technique is particularly costly since it involves a combinatorial optimization over V N , which 
is unfeasible in practice: here the posterior is defined on R 3A and V is the group formed by the permutations 
of K elements, V v has cardinality (Kl) N . As a consequence, in this work, we use the online version of this 
algorithm proposed in 1431 and having a final cost of N(K\). To avoid the use of the resampling in order to 
get this set of unweighted particles, Q, we propose an adaptation of this algorithm, described in Algo. [4] in 
order to be able to use directly the set of weighted particles provided by the SMC sampler. 

IV. Posterior Cramer-Rao bound for multiple source localization 

In this section, we derive the posterior Cramer-Rao bound (PCRB) as an estimation benchmark for the 
parameters. We will thus assume in this setting that we condition on the number of sources. This PCRB 
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Algorithm 4 Online post-processing relabeling algorithm 

1: Input: Set of weighted particles from the last iteration of the SMC sampler Wk*,T 

2: Sort the particle by descending order of their associated weights 
3: Set /ri = e'l'} T and 0^ label = 0^} T and initialize 
4: for n = 2,..., N do 
5: Find 

tin = argmaxA/'(i?(0[” ) T )|At„_i,S„_i') 

V v ’ ' 


6- Set 9^ — ?9 (9^ ) 
O. v^ei V relabel — Vn\V k *,T) 


7: Set = W$ T 


V n 

2^j =1 vv k*,T 


-\ -i 


for i = 1,..., n and compute: 


Pn — U U relabel 


i =1 


- /0(<£Lei - ^) T 


i=l 


8: end for 

9: Use the relabeled collection of weighted particles j 0 


, s —, >-,tv 

( m ) ii/l" 
n vy k 


relabel ’ 
N 


s >1 iV 

% > to compute point estimate, e.g. MMSE: 

’ J m=l 




will thus provide a theoretical performance limit for the Bayesian estimator of the locations as well as the 
transmitted powers of the K sources given the observations, z, obtained at the fusion center. Let us remark 
that in ll25l . the authors have derived the Cramer-Rao bound for the single source problem with quantized 
data and imperfect channel between the sensors and the fusion center. Here, we propose to generalize this 
result by considering Ok as a random variable (Bayesian framework which leads to the posterior CRB) and 
Ok composed of K multiple sources. 

Indeed, the PCRB gives a lower bound for the error covariance matrix Ii44l : 

Tl 


E 


( &k(z ) - Qk) (0k( 


z) - 0 K 


> J 


-l 


(42) 


where J is the 3 K x 3A' Fisher information matrix (FIM) 


J = ^\S7 g ^ogp{z,0K\MK)^^ogp{z,0 K \MK)\ 


= -E 


A® log p(z,O k \Mk) 


(43) 


where Aq := VgV^ is the second derivative operator and V# is the gradient operator with respect to 0. 
Using the fact that p(z, Qk\M-k) = p(z\Qk,-M-k)p(@k\-M-k), the expression of the FIM in ( d43ll ) can be 
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expressed as: 


J = - E 


l\ e e \ogpjz\0 K ,M K ) 


-E 


A ^ \ogp{0 K \M K ] 


— JcL “I - Jpi 


(44) 


where J p represents the a priori information and J,/ is the “standard” FIM (used in the derivation of the 
CRB) averaged over the prior of the different location and power of the K sources: 


Jd= [ Jd{OK)p(0K\MK)dO K - 
Je k 


As demonstrated in the Appendix, this standard FIM is defined for this problem as follows: 


N L—l 


J d ( 0 K ) = 

i =1 3=0 


v ep{zi = j\0 K ,M K )VgP(zi =j\0 K ,M K ) 


p(zi = j\0 K ,M K ) 

with the gradient operator given by: 

V 0 = 

Using ([7]). the gradient term in (l46l) is expressed as: 

L—l 

v e p(zi = j\0 K ,M K ) = £>(* = j\ bi = l )^epjh = l\0 K , M k ), 


8 8 8 

8 Pi dxi dy\ 


8 8 8 
8Pk 8x k dy K \ 


1=0 


in which for k = 1,..., K produces: 

dpjbj = 1\0 k ,M k ) 
dP k 

dpjbj = 1\0 k,M k ) 
dx k 

dpjbj = 1\0k,Mk) 
dy k 


j \ n/2 

«0 \ Pi,l 


do_' n/ 
dj,k, 


dj,k 


and 


Pi,i = I e - e" 


2^27 T0 2 P k 


nPl /2 o 

\ k Pi,l ( Px,i 

- x k ) 


2\/2ito 2 


nPl /2 o 

k,kPhi(Py>i 

- 1Ik) 


2V2iro 2 


( X i,l + l~ a i 

2a 2 

A. 



(45) 


(46) 


(47) 


(48) 


(49) 


(50) 


Although an analytical expression for J,j(O k) has been derived, in order to obtain .J,i involved in the 
computation of the FIM defined in (l44l) . we need to resort to some numerical techniques for the approximation 
of the integral that defines this quantity in (f45l) . The procedure we use is a simple Monte-Carlo integration: 

1) Draw Nmc realization of the state from the prior: ~ p(0/ x ). 

2) Approximate the quantity of interest by: 


-i -* • ivi o 

i=l 


(51) 
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Finally, the second term representing the a priori information in (1441) is a 3K x 2>K matrix defined as: 


a 



0 


with £ = A«+U(o+-i) - Proof: See the Appendix. 



(52) 


V. Numerical Simulations 

In all the experiments, we consider a signal decay exponent and a reference distance as n = 2 and do = 1 
respectively. The ROI is a 100m x 100m. field in which 100 sensors arc deployed in a grid where the location 
of each sensor is assumed to be known. The thresholds of the M-bit quantizer defined in © are the same for 
each sensor and are equally spaced between A,;.i = 0 and Aqi_i = 22, Vi = 1,..., 100. The hyper-parameters 
in the prior distribution of the transmitted power of each source defined in ( flOl ) are a = 50 and b = 2.5 x 10 5 . 
An uniform distribution is used as the prior over the collection of models, i.e. Vfc £ {1,..., K m:ix } we have 
p(Mk) = 1/Amax- All the results have been obtained by using IVmcmc = 5 in the MWG (summarized in 
Algo. El) used in the SMC sampler as forward kernel. In order to illustrate the benefit of using the SMC 
sampler, we have adapted to the problem considered in this paper the importance sampler (IS) that has been 
proposed for a single source localization in lf26l . For each model Aik, this IS algorithm simply consists 
in sampling Njs particles from the prior distribution in ( fTOl ) and in assigning to each of these particles an 
weights which is proportional to the likelihood given in ©. In order to have a fair comparison between 
both algorithms since the proposed SMC sampler adapts the number of iterations on-line using the procedure 
described in Section I1T1-R1 1 the number of particles used in the IS algorithm is set to Njs = TN, with 
T the total number of SMC iterations averaged over multiple runs for the configuration under study. As a 
consequence, the complexity of these two algorithms is equivalent since the number of particles generated 
for all the results described below is the same for both schemes. 


A. Accuracy of the estimators 

We first study the robustness and the accuracy of the estimators of the two main quantities of interest: 
model posterior probabilities and the MMSE of the parameters under each model. In order to perform this 
analysis, both schemes have been run 100 times on the same realization of observations from a scenario 
with 4 sources. From the theory, we know that both IS and SMC algorithms provides, whatever the number 
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of particles used, an unbiased estimate of log p(z\Mk) which corresponds as discussed in dTTb to the only 
unknown quantity in the model posterior distribution. However, as depicted in Fig. [3a] the variance of the 
estimator of this quantity obtained from these two algorithms is significantly different. If both algorithms 
perform similarly for one source, the SMC sampler outperforms significantly the IS algorithm as the number 
of sources increases. The same remark holds for the variance of the MMSE estimator shown in Fig. [3b] which 
is quite remarkable since the MMSE with the SMC sampler is only computed with N particles instead of 
NT with the IS. The same remarks hold even for an increasing number of particles as illustrated in Fig. 
[4] To understand these results, we present in Table Q] the effective sample size (ESS), defined in ([2Tb which 
represents the number of particles that will really contribute to the final estimator. We can see from these ESS 
results that the IS algorithm completely collapses when the dimension of the model (i.e. the number of sources 
in the ROI) increases. As an example, for the model .AT 4 , only 3-4 particles over the A 75 = NT = 1700 
will really contribute to the MMSE estimator by having a non-negligible importance weights. Conversely, the 
ESS from the SMC sampler is quite stable with the dimension of the model. All these results clearly illustrate 
the significant gain provided by the proposed SMC sampler and thus the benefit of gradually introducing the 
likelihood information through the successive iterations of the algorithm. 




(a) (b) 

Fig. 3: Comparison of the variance of both the log evidence, logp(z\Mk), estimator in (a) and MMSE 
estimator in (b) provided by the proposed SMC sampler and the IS algorithm - The parameter of the adaptive 
cooling schedule is set as CESS = 0.9 N leading to the following average number of iterations T,vi, = 8 , 
Tm 2 = 12, Tvf 3 = 15, T_A/f 4 = 17 and T_a ^ 5 = 18 for both N = 50 or N = 100 (Number of quantization 
levels: L = 4 and a 2 = 1) 


Let us now illustrate with a 2 targets scenario, the label switching problem discussed in Section Illl-CI and 
the importance of having a relabeling algorithm in order to provide a point estimate. In Fig. [5] we present 
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Mi 

m 2 

m 3 

Mi 

M 5 

IS 

0.0836 

0.0082 

0.0030 

0.0019 

0.0018 

SMC 

0.6180 

0.6290 

0.6297 

0.6276 

0.6319 


TABLE I: Comparison of the average Effective sample size defined in (I2TT) and scaled by the number of 
particles for both SMC with N = 100 and IS algorithms for the different models 



Fig. 4: Evolution of the variance of the estimators of both the log evidence and the MMSE provided by the 
IS and the SMC for model Mi. For the SMC, the results have been obtained with N = 100 particles by 
varying the CESS thus leading to a different average number of iterations, T. (Number of quantization levels: 
L = 4 and <r 2 = 1 ) 


the marginal posterior distribution obtained with the proposed SMC sampler. We can first remark that the 
algorithm is clearly able to capture the multimodality of each marginals. However, if the MMSE is directly 
computed from this approximation, the estimated //-coordinate for both targets will be approximately 50 
instead of 55 and 45. The proposed relabeling algorithm described in Algo. [4] allows to isolate both modes 
by finding the best permutations for all particles. The MMSE estimate by taking the particle system after 
relabeling will therefore provide an accurate point estimate close to the truth. Moreover from the estimates 
of the posterior distribution in this figure, we can remark that the algorithm is able to detect that there are 2 
targets in the scene. 

Let us now illustrate with Fig. [ 6 J the challenging problem of having two targets of interest that are placed 
very close to each other ([50,49] and [50, 51]). We remark that from the model posterior probabilities that the 
algorithm is able to provide the uncertainty arising from this scenario between the model with 1 and 2 targets. 
The ability of the algorithm to distinguish two close targets will clear be a function of many parameters of the 
sensor network: distance between sensors, variance of the observation noise, number of quantization levels 
as well as the quality of the wireless channel between the sensors and the fusion center. The approximation 
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of the posterior marginal distributions obtained from the SMC sampler under model M 2 are both unimodal 
owing to the relatively small distance between the 2 sources. Once again in this case, we can see the benefit 
of using the relabeling algorithm for computing the final MMSE point estimator. 



V 


(a) Posterior marginal before relabeling 



y 

(b) Posterior marginal after relabeling 


1 



Model 


(c) Model Posterior 

Fig. 5: Illustration of the relabeling and the model posterior obtained with the proposed SMC sampler (N = 
200 and CESS = 0.92/) in a scenario with two targets located at [50,45] and [50, 55] (Number of quantization 
levels: L = 8 and a 2 = 0.5 ) 

Next, we compare the performances of the proposed algorithm when 4 sources are in the ROI. In order 
to obtain the following results, 100 realizations of the four sources and associated observations have been 
drawn from the prior and likelihood defined in Section III-AI Table [TT] illustrates the ability of the proposed 
to detect the correct number of targets. The correct number of target ( i.e. Model M/C) is chosen more often 
with the proposed SMC sampler than the IS algorithm over the 100 scenarios. Even if both provides an 
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(a) Posterior marginal before relabeling 


(b) Posterior marginal after relabeling 



Model 


(c) Model Posterior 


Fig. 6: Illustration of the relabeling and the model posterior obtained with the proposed SMC sampler (N = 
200 and CESS = 0.9JV) in a scenario with two close targets located at [50,49] and [50, 51] (Number of 
quantization levels: L = 8 and cr 2 = 0.5 ) 


unbiased estimate of the evidence of each model, the significant lower variance of the estimator provided 
by the SMC sampler, which was illustrated previously in Fig. 0 allow to have a more efficient and accurate 
model decision step. 
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Mi 

m 2 

M 3 

m 4 

M 5 

IS 

0 

0 

13 

85 

2 

SMC 

0 

0 

3 

96 

1 


TABLE II: Comparison of the number of times that each model has been selected from the estimated model 
posterior probabilities given by both the proposed SMC sampler the IS algorithm under 100 realizations 
(N = 100, CESS = 0.91V, L = 4 and a 2 = 1) 


B. Localization Performance and PCRB 

In Fig- [7] the performance of the proposed SMC sampler (and the IS algorithm) in term of the mean squared 
error between point estimate d p of the algorithm and the true location d p of the four sources: 



MSE = trace <{ 

e ( dp-dp) (d p - 

o P ) T 

} (53) 


1 T 


1 T 


x\ yi ■ 

■ x 4 and d p = 

xi yi ■■■ x 4 

2/4 

represents the estimated (by using 


with d p = 

the relabeling algorithm) and the true location of the four targets, respectively. We also plot the associated 
PCRB that we have derived in Section [TV] In order to obtain the results we use 100 realizations (of the 
different source characteristics and associated observations by avoiding the case in which two targets are very 
close). The results depicted in Figures [7] and [ 8 ] clearly demonstrate the good localization performance of the 
proposed algorithm and the significant gain compared to the IS algorithm which completely fails to localize 
four targets. As expected, the accuracy on the localization improves with the increase of either the number of 
sensors or the number of quantization levels as well as with the decrease of the measurement noise variance. 



Number of quantization levels L 



Number of quantization levels L 


(a) a 2 = 1 


(b) o - 2 = 0.05 


Fig. 7: Evolution of the mean squared error for the source locations as a function of the number of quantization 
levels L with two different values of the measurement noise a 2 (N = 100 and CESS = 0.9A r ) 
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(a) CT 2 = 1 (b) CT 2 = 0.05 

Fig. 8: Evolution of the mean squared error for the source locations as a function of the number of sensors 
in the ROI with two different values of the measurement noise a 2 (Number of quantization levels: L = 8 - 
N = 100 and CESS = 0.91V) 


VI. Conclusion 

In this paper, we addressed the problem of localizing an unknown number of energy emitting sources 
in wireless sensor networks with quantized data. We provided a generalization of recent existing works 
considering a single source. We proposed a Bayesian solution for the joint estimation of the unknown number 
of sources as well as their associated parameters which is based on SMC sampler. Then, we derived the 
posterior Cramer-Rao bound for the estimation of the characteristics of these multiple energy emitting sources. 
Numerical simulations clearly illustrated the ability of the proposed SMC sampler to perform this challenging 
joint estimation. The different experiments showed that the proposed scheme based on novel SMC sampler 
improves quite significantly the accuracy of the estimators method that are required for model selection (i.e., 
the number of sources) and the estimation of the source characteristics compared to more classical importance 
sampling method. 


Appendix 

Proof of the Posterior Cramer-Rao bound 

As presented in Section |IV| the Fisher information matrix (FIM) for the posterior Cramer-Rao bound 
(PCRB) can be decomposed as follows: 


J = 



'•s. 


Jd(^K)p{^K\MK)dd K + E 


Jd 


—Ag \ogp(0 K \M K ) 


(54) 
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where A e e := VgVg is the second derivative operator and Vg is the gradient operator with respect to 9. 
In this appendix, we derive respectively Jd{Oi<) and J p . 

The information matrix J c i(9r) is defined as: 


Jd(9i<) = -^z\0 K Aglogp(z\6 K ,MK) 


(55) 


The first derivative of the log likelihood is given by: 


N 


Vg logp(z\e K ,M r) = ^2V0logp(zi\6 K ,M K ) 


2—1 


VgP(Zi\0 K ,M K ) 
p{zi\6 K: MR) 


E 


(56) 


Therefore, the second derivative can be written as: 


A e e log p{z\9 K , Mr) = V e X7 q log p(z\ 6 k , Mr) 

V"' ^9 V ])P(Zi | 9r , Mr ) 
p{Zi\9 K , Mr) 

yep(zi\9 K ,MR)Vgp{zi\9R,MR) 
p(zi\9 K ,M K ) 2 

To obtain (|55T ). we now take the negative expectation of this second derivative with respect to j)(z,\9 r. Mt = 

K): 


Jd{9R) 


( i a \ a f ^O^gP{ z i ~ j\9R, Mr) 

EEK^*,A(ri{- piZi=jlSk>Mk) 


i= 1 3 =0 

+ 


Vepjzi = j\9R,MR)Vlp(zj = j\9 K ,M K ) 
p(zi = J\9 r,Mr ) 2 


N L—l 

EE 


Vep(zi 


i= 1 3 =0 


j\9R,MR)VgP(Zj = J\9 r,Mr) 
p(zi = j\9 K ,M K ) 


-\7 e \7gp(zi = j\9 K ,M K ). 


(58) 


The second term is equal to 0 since: 

N L—l 

y~! y~! Ve^eP(zi = j\9 K ,MR) 

*=1 3=0 


N L—l 

^ V gVg ^2 p( z i = j\ 9r,Mr) 
i =1 3=0 

" -v-' 

=1 


0. 


(59) 
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As a consequence, we finally obtain: 

T m \ V- V'' v eP(zi = j\e K ,M K )V0p(zi = j\0 K ,M K ) 

‘ ’ = hh --• 

Using (Q, the gradient term involved in this expression can be expressed as: 

L—l 

v 0p{zi = j\0 K , M k ) = £>(* = j\bi = l)V e p(bi = 1\0 k ,Mk), 


(60) 


1=0 


with 


p(bi = i\0 k ,Mk) = Q ( Xhl - Q (Xi ’ l+1 ai 


a 


a 


(61) 


(62) 


As a consequence, since the Q- function is the complementary Gaussian cumulative distribution, we can 
easily remark that: 


= 1\0k,Mk) = 


1 


y/2iro' 2 


e ^ 




VgOj 


(63) 


pi,i 


Finally from the definition of ai in ©, we obtain, for k = I, K: 


dp(bi = 1\6 k ,Mk) 

dP k 

dp(bj = 1\6 k ,M k ) 
dx k 

dp(bj = 1\6 k ,Mk) 

dyk 


do 


n/2 


Pi,l 


di,k / 2\J2 tT(T 2 Pfc 

d 0 \ n / 2 nP\A 2 d~YiAPx,i - x k ) 


i,k 


2y/2i rcr 2 


(64) 


do_\ n/2 nP k /2(i i,kPiAPy,i - Vk) 
di^k, 


2y/2ircr' 2 


which completes the analytical calculation of J ( i{0k)- 

Finally, we derive the a priori information matrix given by: 


J p = E 


-Ag log p(6k\Mk) 


(65) 


From the prior distributions in (IhUTOl). each target’s location and power are independent and identically 
distributed. J p will be therefore a 3 K x 3 K block diagonal matrix with information associated to the location 
and the power defined respectively as: 


E 


~ A \Z’ V yi T lo &N(i x kiyk] T \P'p, £ p ) 



( 66 ) 
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and 


E 


-A%logig(P k \a,b) 


= E 


d 2 

dPi 


log 


r(o) 


P^- 1 exp 


b_ 

P~k 


= E ^( a + 1 ) lo g( p fc) + ^ 

= E[26P,7 3 -(a + l)P,- 2 ] 

= 26E [P- 3 ] - (a + 1)E [P- 2 ] . 


(67) 


Let us now derive the two moments involved in this expression. We have, for n > 0: 


E 


[?-] = [ 
Jo 


+oo ft 

-n 
k 


p ~ n _I p 

r k j-' ' r ' 

foa r+oo 


r(a) 


—a— 1 
k 


exp ( - — ) dP k 


r(a) Jo 

b a r(a + n) 


exp [ ~)dP k 


( 68 ) 


T(a) b a+n 

The last expression is obtained from the expression of the normalizing constant of an inverse-gamma distri¬ 
bution, XQ(a + n,b ). By using the equality of the Gamma function, T(a + 1) = aT(a), we obtain: 


e[p . 2] = (s+i)., 


E[p t : 3 ] = 


b 2 

(cl T- 2)(n *T l)n 

P ' 


(69) 

(70) 


By plugging these expressions in (l67l) . the prior information for the power is given by: 

cl{cl T- l)(n T - 3) 


E 


-A p k k \oglQ(P k \a,b) 


b 2 


= e, 


(71) 


leading to 


Jp — 




Sp 1 . 


(72) 
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